%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cfqdmv}
%
\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pour les notations et l'algorithme dans son ensemble,
on se reportera à \fort{cfbase}.

Dans le premier pas fractionnaire (\fort{cfmsvl}), on a résolu une
équation sur la masse volumique, obtenu une prédiction de la pression
et un flux convectif "acoustique".
On considère ici un second pas fractionnaire au cours duquel seul varie
le vecteur flux de masse $\vect{Q}=\rho\vect{u}$
(seule varie la vitesse au centre des cellules).
On résout l'équation de Navier-Stokes indépendamment
pour chaque direction d'espace, et l'on utilise le flux de masse acoustique
calculé précédemment comme flux convecteur (on pourrait aussi utiliser
le vecteur quantité de mouvement du pas de temps précédent).
De plus, on résout en variable $\vect{u}$ et non $\vect{Q}$.

Le système à résoudre entre $t^*$ et $t^{**}$ est (on exclut
la turbulence, dont le traitement n'a rien de particulier dans le
module compressible)~:

\begin{equation}\label{Cfbl_Cfqdmv_eq_qdm_cfqdmv}
\left\{\begin{array}{l}

\rho^{**}=\rho^{*}=\rho^{n+1}\\
\\
\displaystyle\frac{\partial \rho\vect{u}}{\partial t}+
\divv(\vect{u} \otimes \vect{Q}_{ac}) + \gradv{P}
= \rho \vect{f}_v + \divv(\tens{\Sigma}^v)\\
\\
e^{**}=e^{*}=e^n\\

\end{array}\right.
\end{equation}

La résolution de cette étape est similaire à l'étape
de prédiction des vitesses du schéma de base de \CS.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discrétisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------
\subsection*{Discrétisation en temps}
%---------------------------------

On implicite le terme de convection, éventuellement
le gradient de pression (suivant la valeur de \var{IGRDPP}, en utilisant la
pression prédite lors de l'étape acoustique) et le terme en gradient
du tenseur des contraintes visqueuses.
On explicite les autres termes du tenseur des contraintes visqueuses.
On implicite les forces
volumiques en utilisant $\rho^{n+1}$.

On obtient alors l'équation discrète suivante~:
\begin{equation}\label{Cfbl_Cfqdmv_eq_vitesse_cfqdmv}
\begin{array}{l}
\displaystyle\frac{(\rho\vect{u})^{n+1}-(\rho\vect{u})^n}{\Delta t^n}
+ \divv(\vect{u}^{n+1} \otimes \vect{Q}_{ac}^{n+1})
- \divv(\mu^n \gradt{\vect{u}^{n+1}})\\
\\
\text{\ \ \ \ }= \rho^{n+1} \vect{f}_v - \gradv{\widetilde{P}}
+ \divv\left(\mu^n\ ^t\gradt{\vect{u}^n}
+ (\kappa^n-\frac{2}{3}\mu^n)\divs{\vect{u}^n}\ \tens{Id}\right)\\
\end{array}
\end{equation}
avec $\widetilde{P}=P^n\text{ ou }P^{Pred}$ suivant la valeur de \var{IGRDPP}
($P^n$ par défaut).

En pratique, dans \CS, on résout cette équation en faisant apparaître à
gauche l'écart $\vect{u}^{n+1} - \vect{u}^n$. Pour cela, on écrit la
dérivée en temps discrète sous la forme suivante~:

\begin{equation}
\begin{array}{ll}
\displaystyle
\frac{(\rho \vect{u})^{n+1} - (\rho \vect{u})^n}{\Delta t^n}
& =
\displaystyle
\frac{\rho^{n+1}\, \vect{u}^{n+1} - \rho^n\, \vect{u}^n}{\Delta t^n}\\
& =
\displaystyle
\frac{\rho^{n}\, \vect{u}^{n+1} - \rho^n\, \vect{u}^n}{\Delta t^n}+
\frac{\rho^{n+1}\, \vect{u}^{n+1} - \rho^n\, \vect{u}^{n+1}}{\Delta t^n}\\
& =
\displaystyle
\frac{\rho^{n}}{\Delta t^n}\left(\vect{u}^{n+1} - \vect{u}^n\right)+
\vect{u}^{n+1}\frac{\rho^{n+1} - \rho^n}{\Delta t^n}
\end{array}
\end{equation}

et l'on utilise alors l'équation de la masse discrète pour écrire~:
\begin{equation}
\displaystyle
\frac{(\rho \vect{u})^{n+1} - (\rho \vect{u})^n}{\Delta t^n}
=
\frac{\rho^{n}}{\Delta t^n}\left(\vect{u}^{n+1} - \vect{u}^n\right)-
\vect{u}^{n+1}\dive\,\vect{Q}_{ac}^{n+1}
\end{equation}



%---------------------------------
\subsection*{Discrétisation en espace}
%---------------------------------

%.................................
\subsubsection*{Introduction}
%.................................

On intègre l'équation (\ref{Cfbl_Cfqdmv_eq_vitesse_cfqdmv})
sur la cellule $i$ de volume $\Omega_i$ et
on obtient l'équation discrétisée en espace~:

\begin{equation}\label{Cfbl_Cfqdmv_eq_vitesse_discrete_cfqdmv}
\begin{array}{l}
\displaystyle\frac{\Omega_i}{\Delta t^n}
(\rho_i^{n+1}\vect{u}_i^{n+1}-\rho_i^n\vect{u}_i^n)
+ \displaystyle\sum\limits_{j\in V(i)}
(\vect{u}^{n+1} \otimes \vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}
- \displaystyle\sum\limits_{j\in V(i)}
\left(\mu^n\gradt{\vect{u}^{n+1}}\right)_{ij} \cdot \vect{S}_{ij}\\
\\
= \Omega_i\rho_i^{n+1} {\vect{f}_v}_i
- \Omega_i(\gradv{\widetilde{P}})_i
+ \displaystyle\sum\limits_{j\in V(i)}
\left(\mu^n\ ^t\gradt{\vect{u}^n} + (\kappa^n-\frac{2}{3}\mu^n)
\divs{\vect{u}^n}\ \tens{Id}\right)_{ij}\vect{S}_{ij}\\
\end{array}
\end{equation}

%.................................
\subsubsection*{Discrétisation de la partie ``convective''}
%.................................

La valeur à la face s'écrit~:
\begin{equation}
(\vect{u}^{n+1} \otimes \vect{Q}_{ac}^{n+1})_{ij}\cdot \vect{S}_{ij}
= \vect{u}_{ij}^{n+1}(\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}
\end{equation}
avec un décentrement sur la valeur de $\vect{u}^{n+1}$ aux faces~:
\begin{equation}
\begin{array}{lllll}
\vect{u}_{ij}^{n+1}
& = & \vect{u}_i^{n+1}
& \text{si\ } & (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
& = & \vect{u}_j^{n+1}
& \text{si\ } & (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}
\end{equation}
que l'on peut noter~:
\begin{equation}
\vect{u}_{ij}^{n+1}
 = \beta_{ij}\vect{u}_i^{n+1} + (1-\beta_{ij})\vect{u}_j^{n+1}
\end{equation}
avec
\begin{equation}
\left\{\begin{array}{lll}
\beta_{ij} = 1 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
\beta_{ij} = 0 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}\right.
\end{equation}

%.................................
\subsubsection*{Discrétisation de la partie ``diffusive''}
%.................................

La valeur à la face s'écrit~:
\begin{equation}
\left(\mu^n\gradt{\vect{u}^{n+1}}\right)_{ij}\vect{S}_{ij}
= \mu_{ij}^n
\displaystyle \left( \frac{\partial \vect{u}}{\partial n} \right)^{n+1}_{ij}
S_{ij}
\end{equation}
avec une interpolation linéaire pour $\mu^n$ aux faces (en pratique avec
$\alpha_{ij}=\frac{1}{2}$)~:
\begin{equation}
\mu_{ij}^n
= \alpha_{ij}\mu_{i}^n+(1-\alpha_{ij})\mu_{j}^n
\end{equation}
et un schéma centré pour le gradient normal aux faces~:
\begin{equation}
\displaystyle \left( \frac{\partial \vect{u}}{\partial n} \right)^{n+1}_{ij}
= \displaystyle\frac{\vect{u}_{J'}^{n+1}-\vect{u}_{I'}^{n+1}}{\overline{I'J'}}
\end{equation}

%.................................
\subsubsection*{Discrétisation du gradient de pression}
%.................................

On utilise \fort{grdcel} standard. Suivant la valeur de \var{IMRGRA},
cela correspond à une reconstruction itérative ou par moindres carrés.

%.................................
\subsubsection*{Discrétisation du ``reste'' du tenseur des contraintes visqueuses}
%.................................

On calcule des gradients aux cellules et on utilise une
interpolation linéaire aux
faces (avec, en pratique, $\alpha_{ij}=\frac{1}{2}$)~:
\begin{equation}
\begin{array}{r}
\left(\mu^n\ ^t\gradt{\vect{u}^n} + (\kappa^n-\frac{2}{3}\mu^n)
\divs{\vect{u}^n}\ \tens{Id}\right)_{ij}\cdot\vect{S}_{ij}
= \left\{\alpha_{ij} \left(\mu^n\ ^t\gradt{\vect{u}^n}
+ (\kappa^n-\frac{2}{3}\mu^n)\divs{\vect{u}^n}\ \tens{Id}\right)_i\right.\\
\\
\left.+ (1-\alpha_{ij}) \left(\mu^n\ ^t\gradt{\vect{u}^n}
+ (\kappa^n-\frac{2}{3}\mu^n)\divs{\vect{u}^n}\ \tens{Id}\right)_j
\right\} \cdot\vect{S}_{ij}\\
\end{array}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On résout les trois directions d'espace du système
(\ref{Cfbl_Cfqdmv_eq_vitesse_discrete_cfqdmv}) successivement et indépendamment~:
\begin{equation}\label{Cfbl_Cfqdmv_eq_vitesse_finale_cfqdmv}
\left\{\begin{array}{l}
\displaystyle\frac{\Omega_i}{\Delta t^n}
(\rho_i^{n+1}{u_i}_{(\alpha)}^{n+1}-\rho_i^n{u_i}_{(\alpha)}^n)
+ \displaystyle\sum\limits_{j\in V(i)}
{u_{ij}}_{(\alpha)}^{n+1}(\vect{Q}_{ac}^{n+1})_{ij}\cdot\vect{S}_{ij}
- \displaystyle\sum\limits_{j\in V(i)}
\mu_{ij}^n\frac{{u_j}_{(\alpha)}^{n+1}-{u_i}_{(\alpha)}^{n+1}}{\overline{I'J'}}S_{ij}\\
\qquad\qquad\qquad\qquad= \Omega_i\rho_i^{n+1} {{f_v}_i}_{(\alpha)}
- \Omega_i{(\gradv{\widetilde{P}})_{i}}_{(\alpha)}\\
\qquad\qquad\qquad\qquad + \displaystyle\sum\limits_{j\in V(i)}
\left((\mu^n\ ^t\gradt{\vect{u}^n})_{ij}\cdot\vect{S}_{ij}\right)_{(\alpha)}
 + \displaystyle\sum\limits_{j\in V(i)} \left((\kappa^n-\frac{2}{3}\mu^n)
\divs{\vect{u}^n}\right)_{ij}{S_{ij}}_{(\alpha)}\\
i = 1\ldots N \qquad \text{et} \qquad (\alpha) = x,\ y,\ z\\
\end{array}\right.
\end{equation}

Chaque système associé à une direction est résolu par une méthode
d'incrément et résidu en utilisant une méthode de Jacobi.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Points à traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% propose en patch 1.2.1

%Compléter le commentaire en entête de \fort{visecv} pour prendre en compte
%la viscosité en volume.
